library(tidyverse)
library(sf)
library(mapview)
library(leafpop)
library(tigris)
library(DT)
starbucksNC <- read_csv("https://raw.githubusercontent.com/libjohn/mapping-with-R/master/data/All_Starbucks_Locations_in_the_US_-_Map.csv") %>%
filter(State == "NC")
## Parsed with column specification:
## cols(
## .default = col_character(),
## `Store Number` = col_double(),
## `Facility ID` = col_double(),
## `Food Region` = col_double(),
## Latitude = col_double(),
## Longitude = col_double()
## )
## See spec(...) for full column specifications.
options(tigris_use_cache = TRUE)
durham_blocks <- tracts(state = "NC", county = "Durham", class = "sf", year = 2012)
durham_blocks <- durham_blocks %>%
select(-6, -7, -8, -11, -12)
glimpse(durham_blocks)
## Observations: 60
## Variables: 8
## $ STATEFP <chr> "37", "37", "37", "37", "37", "37", "37", "37", "37",...
## $ COUNTYFP <chr> "063", "063", "063", "063", "063", "063", "063", "063...
## $ TRACTCE <chr> "001001", "000200", "000900", "000302", "000102", "00...
## $ GEOID <chr> "37063001001", "37063000200", "37063000900", "3706300...
## $ NAME <chr> "10.01", "2", "9", "3.02", "1.02", "10.02", "4.01", "...
## $ ALAND <dbl> 2283723, 3021358, 1025539, 1722150, 4011880, 3348806,...
## $ AWATER <dbl> 3268, 0, 0, 0, 3346, 0, 43208, 0, 121500, 0, 3875, 18...
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-78.88456 3..., MULTIPOL...
sbux_drm <- starbucksNC %>%
filter(City == "Durham") %>%
select(`Store Number`, `Ownership Type`,
`Phone Number`, City, Zip, Longitude,
Latitude)
# assign my coordinates to a degree system
# documenting my XY values in degrees
# use, for example 4326 (best for US) ; 4269 (for US Agency)
sbux_drm_sf <- st_as_sf(sbux_drm,
coords = c("Longitude", "Latitude"),
crs = 4326)
st_crs(starbucksNC)
## Coordinate Reference System: NA
st_crs(sbux_drm_sf)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(durham_blocks)
## Coordinate Reference System:
## EPSG: 4269
## proj4string: "+proj=longlat +datum=NAD83 +no_defs"
# https://guides.library.duke.edu/r-geospatial/CRS
# > Projections / CRS > ESPG codes for CR's
Use one of the NC Projections
# Transform to a Projection
durham_blocks <- st_transform(durham_blocks,
st_crs(2264))
sbux_drm_sf <- st_transform(sbux_drm_sf,
st_crs(durham_blocks))
st_crs(sbux_drm_sf)
## Coordinate Reference System:
## EPSG: 2264
## proj4string: "+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"
st_crs(durham_blocks)
## Coordinate Reference System:
## EPSG: 2264
## proj4string: "+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"
mapview(durham_blocks)
mapview(sbux_drm_sf)
mapview(list(durham_blocks, sbux_drm_sf))
shp_with_bux <- st_join(durham_blocks, sbux_drm_sf, join = st_contains, left = FALSE)
bux_within_shp <- st_join(sbux_drm_sf, durham_blocks, left = TRUE, join = st_within)
plot(st_geometry(shp_with_bux))
plot(st_geometry(bux_within_shp))
mapview(durham_blocks, col.regions = "grey",
legend = FALSE,
popup = popupTable(st_drop_geometry(durham_blocks),
zcol = c(4,5),
row.numbers = FALSE,
feature.id = FALSE)) +
mapview(shp_with_bux,
zcol = "NAME",
layer.name = "Tract Name") +
mapview(bux_within_shp,
legend = FALSE,
cex = 0,
#cex = "AWATER",
lwd = 3,
alpha = .5,
color = "red",
zcol = "Store Number",
popup =
popupTable(st_drop_geometry(bux_within_shp),
zcol = c(1,2,3,5,10,11),
row.numbers = FALSE,
feature.id = FALSE))
datatable(st_drop_geometry(bux_within_shp))
bux_dt <- st_drop_geometry(bux_within_shp) %>%
select(1, GEOID, TractName = NAME, Ownership = `Ownership Type`) %>%
group_by(TractName) %>%
summarise(TotalStores = n()) %>%
arrange(-TotalStores)
datatable(bux_dt)
Shade Tract shape by number of Stores in the Tract
bux_tot_in_shape <- left_join(shp_with_bux, bux_dt, by = c("NAME" = "TractName"))
mapview(durham_blocks, col.regions = "grey",
legend = FALSE,
popup = popupTable(st_drop_geometry(durham_blocks),
zcol = c(4,5),
row.numbers = FALSE,
feature.id = FALSE)) +
mapview(bux_tot_in_shape,
zcol = "TotalStores",
layer.name = "Total Stores by Tract",
col.regions = c("royalblue1", "royalblue4")) +
mapview(bux_within_shp,
legend = FALSE,
cex = 0,
#cex = "AWATER",
lwd = 4,
color = "yellow",
zcol = "Store Number",
popup =
popupTable(st_drop_geometry(bux_within_shp),
zcol = c(1,2,3,5,10,11),
row.numbers = FALSE,
feature.id = FALSE))